We chose α-Fe with a body-centered cubic (bcc) structure to develop the coarse-grained model of Fe. Figure 2 depicts the mapping scheme of the α-Fe from the all-atom model to bead–spring model. Figure 2(a) shows an all-atom model of the α-Fe built using Materials Studio (MS) software. The dimension of the simulation box is 22.93 ų. It contains 512 crystal lattices with a lattice parameter of 2.86 Å. We applied the periodic boundary conditions to the simulation box. The COMPASS (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies) forcefield and the Nosé–Hoover thermostat (∼ 300 K) were employed during simulation. The system reached the equilibrium state after running for 50 ps with a time step of 1 fs, as shown in Fig. 2(a). When we built the coarse-grained model of the α-Fe, one α-Fe crystal lattice was treated as one Fe bead, and its geometric center was defined as the coordinate of the corresponding bead. Then, the radial distribution functions (RDFs) of the Fe bead could be calculated on the basis of the equilibrium distribution of the α-Fe atoms, as shown in Fig. 2(b). The RDF of Fe was considered as the target distribution of the coarse-grained Fe model. Figure 2(c) shows the initial regular distribution of the 512 Fe beads built in the large-scale atomic/molecular massively parallel simulator (LAMMPS).^[26] The dimensions of the simulation box are 22.93 Å×22.93 Å×21.06 Å. A Lennard–Jones truncated and shifted (LJTS) potential was used to describe the interaction between the Fe beads; its mathematical expression is as follows:^[27]

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Mapping scheme of the α-Fe from the all-atom model to coarse-grained model: (a) all-atom model of α-Fe atoms, (b) comparison of the RDF distributions, (c) regular distribution, (d) random distribution, and (e) equilibrium distribution of the Fe beads.

To run the coarse-grained MD simulation, firstly, we initialized ε and σ. Based on the mapping scheme of α-Fe from the all-atom model to bead spring model, the central distance between two neighboring beads (2.8 Å) was set to be the initial value of σ. Considering that the interaction strength between the Fe beads is greater than that between the PTFE beads,^[23] the initial value of ε was chosen to be 10 kcal/mol. Then, the size of the coarse-grained model was adjusted using the “fix deform_box” command^[26] until it was similar to that of the all-atom model. The simulation was conducted in the Langevin thermostat^[28] assuming a microcanonical ensemble (also called NVE ensemble). The size of the system reached 22.93 ų after running for 10 ps with a time step of 1 fs, as shown in Fig. 2(d). The coarse-grained MD simulation acquired an equilibrium state after running for 50 ps by maintaining a constant box size, Langevin thermostat, and time step, as shown in Fig. 2(e).

The RDF distribution g_i(r) obtained from the coarse-grained Fe model based on the estimated ε and σ might be different from the target RDF distribution g_target(r). Therefore, we chose to employ iterative Boltzmann inversion (IBI) method^[29] to optimize and determine the final values of ε and σ by comparing the difference in RDFs acquired from the all-atom model and coarse-grained model, respectively. This process can be expressed as follows:^[29]

The values of ε and σ for the coarse-grained Fe model were determined to be 80 kcal/mol and 2.45 Å, respectively, after a number of iterations. A comparison between the RDF distributions of the coarse-grained model and the all-atom model is shown in Fig. 2(b). We can observe that both the RDF distributions of the coarse-grained model and all-atom model contain a number of peaks, and the positions of their first peaks are well in line with each other. Although a slight RDF distribution offset is observed from the second peak onward, the overall trend of these two distributions remains the same. Considering that the interaction between the two Fe beads sharply decreases with an increase in their separation distance, we chose 80 kcal/mol and 2.45 Å, respectively, as the final values for ε and σ to describe the interaction between the Fe beads in the developed coarse-grained model. While developing the coarse-grained model of PTFE, the PTFE moelcule with a molecular weight of 1038 was buit.^[30] It contained ten repeat units and one repeat unit was simplied into one bead.^[23] Therefore, one PTFE molecular chain was simplied into eight backbone beads and two end beads. The interaction parameters among the PTFE beads were depicted in our previous work^[23] in detail.

The interactions between Fe beads and PTFE beads are described while simulating the friction process between Fe and PTFE. In this study, the LJ potential was adopted to describe the interactions between the Fe and PTFE beads, and the combination rule was used to determine the Fe–PTFE interaction parameters as follows:^[31]

When simulating the Fe–PTFE friction process, we defined three kinds of non-bonding interactions; the interaction between the Fe beads was called non-bonding_1, the interactions among the PTFE beads were called non-bonding_2, and the interactions between the Fe beads and PTFE beads were called non-bonding_3, as shown in Fig. 3.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Schematic of the interactions between beads.

To simulate the friction process between Fe and PTFE, we built a two-layer friction model including Fe and PTFE layers. Figure 4 depicts the detailed modeling process. Figure 4(a) shows the regular initial distribution of the Fe and PTFE beads. The upper Fe layer with a size of 30 Å×32 Å×240 Å contains 8000 Fe beads. The lower PTFE layer with a size of 216 Å×32 Å×276 Å contains 1600 PTFE molecules. Considering that the size of the upper specimen is usually smaller than that of the bottom one in the contact surface during a friction test, the length of the upper Fe layer is smaller than that of the bottom PTFE layer along the x direction in our MD simulation. The relaxation process of the Fe–PTFE friction model is divided into two steps.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Fe–PTFE friction model: (a) initial distribution of the Fe beads and PTFE beads, (b) random distribution of the PTFE beads, (c) compressed PTFE molecules, (d) equilibrium distribution of the PTFE molecules, (e) compressed Fe layer, (f) equilibrium distribution of the Fe beads, (g) final distribution of the Fe and PTFE beads.

In the first step, we ran the PTFE layer until it reached an equilibrium state. To decrease the effect of the regular initial configuration on the simulation results, we adjusted the cutoff distance of the LJ potential (non-bonding_2) until only repulsive forces existed between the PTFE beads. However, the distribution controlled by the repulsive forces cannot be random for a short running time. When run for a long time, spacing between the PTFE molecules is large and they scatter in the simulation box. To randomize the PTFE molecules within a certain space, the “fix wall/reflect” command was used to restrict the motion of the PTFE molecules in the range of 0–280 Å along the z direction. Periodic boundary conditions were applied in x and y directions to remove edge effects. The simulation was performed assuming a canonical ensemble (NVT ensemble) with a time step of 1 fs. The PTFE molecules were randomly distributed after running for 50 ps, as shown in Fig. 4(b). To bring the density of the PTFE layer closer to reality (∼ 2.2 g/cm³), we compressed the PTFE layer along the z direction for 30 ps (Fig. 4(c)). Then, the simulation was continued for 20 ps until the PTFE layer reached equilibrium. The Langevin thermostat considers the effect of friction force between the PTFE beads when controlling the system temperature. Thus, the Langevin thermostat was used to control the temperature based on the microcanonical ensemble (NVE ensemble) instead of the NVT ensemble. The equilibrium distribution of the PTFE molecules was obtained after 50 ps, as shown in Fig. 4(d). The dimension of the PTFE layer along the z direction is approximately 185 Å.

In the second step, the relaxation of the Fe layer was conducted. Considering that the simulation would be interrupted if the initial distance between the Fe beads was too short, we assigned a relatively large (∼ 3 Å) initial spacing between the Fe beads, and then slowly compressed the Fe beads along the x direction. The non-bonding interaction (non-bonding_1) would occur in the Fe beads while maintaining a certain distance. Figure 4(e) shows the distribution of the Fe beads after being compressed for 20 ps. We can observe that the Fe layer deforms and shrinks in the x and z directions from 30 Å and 240 Å to 22 Å and 231 Å, respectively, due to the non-bonding interaction between the Fe beads. Based on the Langevin thermostat and NVE ensemble, the relaxation was sustained for 30 ps. The equilibrium distribution of the Fe beads is shown in Fig. 4(f).

To enable the friction model to reach equilibrium, we continued to run the simulation for 50 ps. The final distributions of the Fe and PTFE beads are shown in Fig. 4(g). The bottom of the PTFE layer (z < 20 Å) was fixed to provide support to the system during simulation. Moreover, to reduce the effect of the initial contact position on the final simulation results, three contact positions (Fig. 4(g)) were randomly selected. Further, the average of the simulation results for the selected contact positions was considered as the result of the corresponding condition.

Figure 5 shows the upper Fe layer in contact with the bottom PTFE layer at position two. When simulating the friction process, we first reduce the spacing between Fe and PTFE to approximately 5 Å, as shown in Fig. 5(a). Then, the Fe layer with an external load of 9 MPa is moved along the x direction at a velocity of 0.05 m/s. The Fe layer starts to come in contact with and immerses into the PTFE layer, and the immersion depth increases with an increase in the simulation time. The distributions of the Fe and PTFE layers, after 200 ps, are shown in Fig. 5(b).

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Contact and separation processes between Fe and PTFE layers when the external load is 9 MPa and the velocity is 0.05 m/s: (a) T = 0 ps, (b) T = 200 ps.

By varying the external load applied on the Fe layer, we can obtain the effect of the external load on the friction coefficient of Fe–PTFE at a fixed velocity of 0.05 m/s. Figure 6 shows a comparison between the MD simulation and experimental results. We can see that the friction coefficient acquired from the MD simulation decreases linearly with the external load increasing, which is in good agreement with the experimental results. To elucidate the mechanisms of influence of the external load on the friction coefficient of PTFE, we analyze the variations of radius of gyration R_g, energies of the PTFE layer, and the interaction forces between the Fe and PTFE layers with different external loads.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Comparison between experimental and MD simulation results.

The deformation of polymers can be described using radius of gyration R_g. A high probability density of R_g indicates that more polymer molecules exist in the corresponding shape. Here, we calculate the x and z components of R_g on the basis of the following equation^[32] to elucidate the deformation of PTFE molecules under different external loads:

The distributions of the x and z components of R_g are plotted in Fig. 7. Here, we consider the distributions of R_gx and R_gz corresponding to the absence of contact between the Fe to PTFE as reference lines. We note that the peaks of R_gx and R_gz are in the ranges of 6–8 Å and 1–3 Å, respectively, by comparing Fig. 7(a) and Fig. 7(b). This indicates that most of the PTFE molecules are distributed along the x direction and slightly bended along the z direction. With an increase in the external load, the probability densities of R_gx and R_gz located in the ranges of 6–8 Å and 1–3 Å, respectively, reduce. This implies that a portion of the PTFE molecules experience severe bending along the x and z directions, and the degree of bending rises with an increase in the external load. Moreover, by comparing Fig. 7(a) and Fig. 7(b), we can observe that more severe PTFE molecules deformation occurs in the x direction than that in the z direction. This is probably caused by the relatively low stiffness of the PTFE molecules along the main-chain direction. Figure 8 depicts the structural evolution of one PTFE molecular chain for different loads during the friction processes. We can observe that the initial distribution the PTFE molecule is along the x direction and slightly bends along the z direction, which corresponds to the results shown in Fig. 7. With an increase in simulation time and external load, a distinct bend is observed. This is consistent with the variation trends of R_gx and R_gz (Fig. 7).

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Distributions of (a) the x and (b) z components of the radius of gyration.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Snapshots of the structural evolution of one PTFE molecular chain under (a) 5 MPa, (b) 11 MPa, and (c) 18 MPa during the friction processes.

Figure 9 shows the following energy evolutions for different external loads: bond energy, bond angle energy, kinetic energy of the PTFE layer, and the interaction energy between the Fe and PTFE layers. When the Fe layer starts to come in contact with the PTFE layer, no variation can be observed in terms of the energy distributions for different external loads. However, with an increase in the contact time, the bond energy and bond angle energy gradually increase, as shown in Figs. 9(a) and 9(b). The differences of the bond energy and bond angle energy gradually become larger for different loads. This is because a rise in the external load leads to the bond length and bond angle deviating from their equilibrium values, thereby enhancing the bond energy and bond angle energy. The variations of the bond length and bond angle also result in the deformation of the PTFE molecules (Fig. 8), which corresponds to the changes of R_gx and R_gz (Fig. 7). The kinetic energy of the PTFE layer slightly fluctuates during contact for different external loads. Further, the kinetic energy increases a bit as the external load increases as shown in Fig. 9(c). With an increase in the load, the Fe layer immerses deeper into the PTFE layer within the same contact time, thereby leading to more PTFE molecules moving far away from their original positions, in turn, enhancing the velocity of the PTFE beads, and finally generating a larger kinetic energy. Figure 9(d) shows the non-bonding interaction energy between the Fe and PTFE layers. The minus sign denotes a repulsive force. We see that the non-bonding interaction energy increases with an increase in the contact time and external load.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (a) Bond energy, (b) angle energy, (c) kinetic energy, and (d) Fe–PTFE interaction potential energy as a function of simulation time for different external loads.

We also analyze the effect of the external load on the friction force and normal force, as shown in Fig. 10. It can be observed that both the friction force and normal force increase with the external load increasing. Moreover, the normal force has a higher increment rate than the friction force, consequently reducing the friction coefficient (Fig. 6). During the Fe to PTFE contact, the friction force and normal force are primarily attributed to the variation of the energies shown in Fig. 9. Therefore, by combining Figs. 9 and 10, we arrive at a few potential cause for the change in the Fe–PTFE interaction forces with the external load. Firstly, the higher external load results in a larger immersion depth during the same contact time, thereby increasing the Fe–PTFE contact area, ultimately enhancing the non-bonding interaction force between the Fe and PTFE layers. Secondly, for higher external load, the greater degree of springback caused by the severely deformed PTFE molecules also contributes to an increase in the friction force and normal force. Moreover, the stiffness of the PTFE molecules along the main-chain direction is lower than that along the z direction, resulting in an increment of the friction force that is lower than that of the normal force for different external loads, thus generating a lower friction coefficient for a higher external load. Thirdly, the increase of the kinetic energy also contributes to the Fe to PTFE interaction forces. A larger kinetic energy would increase the Fe–PTFE interaction probability, thereby raising the friction force and normal force.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. Variation of (a) friction force and (b) normal force with simulation time under different loads.